function [S,F] = ident_ext_instr(n_vars,p,instr,res)

% Perform the identification routine
% Inputs are the vector of dependent variables (Y), the number of lags (p),
% the reduced form residuals, the instruments, and the starting dates
% The output is the vector of structural responses S, F-stat for 1st stage,

res_temp = res;

% Compute the Sigma matrix, covariance of reduced form residuals
var_sigma = res_temp'*res_temp/(size(res_temp,1) - n_vars*p - 1);

% 1st stage regression: regress the reduced form residuals on the surprises
Phib  = [ones(size(instr,1),1), instr]\res_temp; %betas
uhat  = [ones(size(instr,1),1), instr]*Phib;     %forecasts, uhats
uhatp = uhat(:,1);                               %the uhat for the policy instrument
rss_1 = sum((res_temp(:,1) - uhat(:,1)).^2);     %RSS for the large model

% Compute the F-statistic for the 1st Stage regression
% Need to compare the SSE for a regression with intercept only
Phib_inter = ones(size(instr,1),1)\res_temp; %intercept only
uhat_0     = ones(size(instr,1),1)*Phib_inter; %fitted values from using intercept only
rss_0      = sum((res_temp(:,1) - uhat_0(:,1)).^2); %RSS for the small model

F_n = (rss_0 - rss_1)/(size(instr,2));              %numerator for F-statistic
F_d = rss_1/(length(res_temp) - 1 - size(instr,2)); %denominator for F-statistic
F   = F_n/F_d; %F-stat for 1st stage regression

% Regress the reduced form residuals on the u_hat_p
gam = [ones(size(uhatp,1),1), uhatp]\res_temp; %this is the sq/sp
gam = gam(2,2:end)';                           %exclude the sp/sp element

% Use the variance-covariance matrix to identify sp, this part follows the
% footnote in pg. 13 of the working paper
% Only difference in notation is s21/s11, which I call "gam"
Sig11 = var_sigma(1,1);
Sig21 = var_sigma(2:end,1);
Sig22 = var_sigma(2:end,2:end);

Q        = gam*Sig11*gam' - (Sig21*gam' + gam*Sig21') + Sig22; %line (18)
s12s12p  = ((Sig21 - gam*Sig11)'/Q)*(Sig21 - gam*Sig11);       %line (17)
s11      = sqrt(Sig11 - s12s12p);                              %line (16)

% Finally, construct the S vector
S = [s11; s11*gam];

end